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We study the spin transport in the quantum Heisenberg spin chain subject to boundary magnetic 
fields and driven out of equilibrium by Lindblad dissipators. An exact solution is given in terms of 
matrix product states, which allows us to calculate exactly the spin current for any chain size. It 
is found that the system undergoes a discontinuous spin-valve-like quantum phase transition from 
ballistic to sub-diffusive spin current, depending on the value of the boundary fields. Thus, the 
chain behaves as an extremely sensitive magnonic logic gate operating with the boundary fields as 
the base element. 


I. INTRODUCTION 

One of the fundamental issues in condensed matter 
physics is the determination of macroscopic parameters 
from the underlying microscopic properties. For systems 
in equilibrium, the Gibbsian approach gives an elegant 
solution since it depends only on the underlying micro¬ 
scopic energy spectrum. However, even if substantial 
progress has recently been made in understanding non¬ 
equilibrium systems, in particular through the so called 
fluctuation theoremSfi^— no such approach is available 
for systems in a Non-Equilibrium Steady-State (NESS), 
characterized by the existence of steady currents. This 
forces one to resort to a full dynamical calculation in 
order to extract steady-state parameters. Such a diffi¬ 
culty is inherent of non-equilibrium systems, dating back 
to Drude’s calculation of the electrical conductivity of 
metals in 1900.— As another example, we note the re¬ 
cent discussions concerning the microscopic derivation of 
Fourier’s law in insulating crystals— “— 

A more thorough understanding of the NESS is also 
essential for the development of several applications in 
phononicSfi^^— spintronicsf^^^— and magnonics— We 
point in particular to two recent remarkable pap ers by 
Ghumak et. al^ and Oltscher et. ah^. In Ref. l20l the au¬ 
thors report on a magnonic logic gate, where the magnon 
current is adjusted by controlling the number of magnon 
scattering processes induced by an auxiliary magnon in- 
jector (the base). On a different setting the authors in 
Ref. study the transport of spin polarized current in 
a two dimensional electron gas. They observe for the 
first time the existence of a ballistic spin flow, in stark 
disagreement with classical predictions. 

The transport properties reported in Refs. [I^ and [13 
both involve the presence of a NESS. Moreover, they 
share in common the fact that they cannot be explained 
by classical theories, thus requiring a full quantum treat¬ 
ment. On the theoretical side, these quantum NESS 
are usually implemented on Id lattice spin systems cou¬ 
pled to external reservoirs j ^ The effect of the reser¬ 
voirs is quite often described by a non-unitary Lindblad 


dynamical equation— However, these models, being 
quantum many-body problems, can seldom be solved ex¬ 
actly and from a numerical point of view they can usually 
only be solved for small lattices. 

The purpose of this paper is to study the transport 
properties in the NESS of the one-dimensional Heisen¬ 
berg chain coupled to two Lindblad reservoirs at each 
end, and also subject to magnetic fields at its bound¬ 
aries. Remarkably, the steady state of this model is ex¬ 
actly expressible in terms of a matrix product stat o^^i^^ 
involving operators satisfying the SU(2) algebra (in the 
case of an XXZ chain this generalizes to the quantum 
U,[SU(2)] algebra31). This provides a method to com¬ 
pute the steady-state spin current J for any chain size^ 
We will show that depending on the strength of the ap¬ 
plied magnetic field, J may undergo a discontinuous spin- 
valve-like quantum phase transition from ballistic to sub- 
diffusive (J ^ 1/AI^; cf. Fig.[3Kd) below). As we shall dis¬ 
cuss, the origin of this transition is related to the entrap¬ 
ment of magnons inside the chain caused by the bound¬ 
ary fields which, in turn, increase the number of magnon 
scattering events. We argue that our system may be used 
as an extremely sensitive magnonic logic gate operating 
with an external magnetic field as the base element. 


II. DESCRIPTION OF THE MODEL 

We consider the isotropic Heisenberg spin-1/2 chain 
with N sites described by the Hamiltonian 

JV-l 

^ = 2 E + «+i) + 

( 1 ) 

where the cr’s are the usual Pauli matrices. The last 
term describes the Zeeman interaction experienced by the 
boundary spins with a field pointing in the 2 ; direction on 
the first site and in the —n„ = (sin 9, 0, — cos 9) direction 
on the last site. Note that with this parametrization the 
boundary fields point in opposite directions when 9 = 0. 
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The chain is coupled to two reservoirs at each end such 
that its density matrix p is governed by the Lindblad 
master equatioit^ 

^ = _,[H,p]+DUp)+Di,ip), (2) 

where the left and right dissipators Dl{r) are given by 
D^{p) = ^ - {K^K,p}, ( 3 ) 

r—zb 

with = 1 / 7(1 ± f)at and =F and 

where are the ladder operators in the direction. 
Explicitly one has = (cosdtr^ — icr^ + sin0cr^)/2 
for the lowering operator and the adjoint expression for 
the raising one. 

The forcing term / S [0,1] describes the polarization 
of the spin reservoirs and is related to a reservoir inverse 
temperature /3 by / = tanh(/l). At / = 1 (zero tem¬ 
perature) the left bath corresponds to a perfect magnon 
source, pumping magnons into the system at a rate 7 , 
while the right dissipator is a perfect drain, absorbing 
magnons at the same rate. We shall concentrate mostly 
on / = 1 , even though some words will be given for / < 1 . 
Note that when h > 0 {< 0) the boundary fields point in 
the same (opposite) direction as the dissipators (Fig.[T|). 
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FIG. 1. (Color Online) Schematic drawing of the dissipators 
Dl,r{p) (black, solid arrows) and the boundary fields h (red, 
dashed arrows) acting on the first and last spins of the chain. 
In (a) the fields are in the same direction as the dissipators 
{h > 0 ) and in (b) the fields act in directions opposite to the 
dissipators (h < 0 ). 


III. OUTLINE OF THE MATRIX PRODUCT 
STATE SOLUTION 

The unique NESS attained by the system at long times 
is the solution of Eq. ([2]) with dp/ dt = 0: 

i[H, p] = Dl{p) + Dr{p) . (4) 

At / = 1 the exact solution was found in Ref. [23 in terms 
of a matrix product state (MPA), as we now outline. The 


first step is to note that since p is a Hermitian positive 
semi-definite operator, we may use the following param¬ 
eterization: 

SS^ 

^ tr(,S,St) ■ ( ) 

For a Heisenberg chain made of N spins 1/2, the operator 
S lives on the Hilbert space Sj = . 

We now use the ansatz that S can be described by a 
matrix-product state: 

s = ( 6 ) 


where H is a 2 x 2 matrix with operator-valued entries 
( Sz S+ \ 

0=( ^ j = -h 5 '+(t+- h S'-CT” . (7) 

The operators Sa live in an auxiliary space 21 so that 
G (g) 21. After contracting with |(/) and \il}) we 
recover 5 G ij. From the bulk structure of the Hamilto¬ 
nian (ED, it can be shown that if Eq. ([HD is to be a solution, 
then the operators Sa must obey the SU(2) algebra: 


[Sz,S±] = ±S± 

( 8 ) 

[S+,S.] = 2Sz. 

(9) 


The proper representation of the algebra to be explic¬ 
itly used in the MPS solution is specified by a complex 
representation parameter p which is fixed by substitut¬ 
ing Eq. dnD and into the steady-state equation o 
and solving the resulting equations. It turns out that 
it is fixed by a lowest weight condition: {4'\Sz = p(</|, 
where {<j)\ = ( 0 | is a lowest weight state of the represen¬ 
tation. Explicitly, in terms of a semi-infinite set of states 
{|n)}()To one has the irreducible representations 


00 


Sz = 

( 10 ) 

n —0 


00 

S+ = ^(n-b l)|n)(n-b 1| 

( 11 ) 

n —0 


00 

S- = '^{2p-n)\n + l){n\ . 

n —0 

( 12 ) 


Notice that for half-integer values of p these representa¬ 
tions reduce to the usual finite dimensional representa¬ 
tions of SU(2). In the present case, the representation 
parameter turns out to be 


i 

2(7 — ih) 


(13) 


which fixes the associated infinite dimensional represen¬ 
tation of SU(2). The right state \ip) over which is 
evaluated is given by the coherent state 


1^) = V’= - tan(6»/2) . (14) 

n=o V” / 
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Including these results in Eqs. m and (ED gives a com¬ 
plete solution for the density matrix of the steady-state. 

From this general solution it is possible to compute 
the expectation value of any local observable;^ the most 
important of which is the spin current Ji leaving site 
i toward site z -I- 1. It is defined from the continuity 
equation 


dt 


— Ji—l Ji 


where 




These equations are valid for * = 2,..., TV — 1. Slightly 
different equations apply to the boundaries. In the steady 
state d(crf)/dt = 0 which gives 


where Z{N) is the normalization constant, 

Z{N) = tT{p) = {0,0\B^\^|;,r) ■ ( 20 ) 

The explicit computation of Z[N) requires constructing 
the matrix 

iBo)k,e = 2|p — k\‘^Sk,e + + |2p — 

We then have 


N 


Z(TV) = ^«)o.fctan 2 '“'(^/ 2 ) 




In particular, when ^ = 0 one has simply Z{N) = 
(So^)o.o. 

It can be further shown that 


Ji = J 2 = ... = Jn '■= J 


[B,,By]=2t{T,-S,)Bo 


The expectation value of an arbitrary observable A 
may be computed as 


(A) = ti{Ap) 


tr(S'US') 

tr(5tS') 


Our strategy will be to first trace over the Hilbert space 
and write everything in terms of expectation values on 
the auxiliary space. But note that S and will each 
contain an auxiliary space. So when we write SS^ we 
must double our auxiliary space. That is, we write 

55t = (0,0|O(p)HT(p*)|V;,zA*) 

where H(p) and act on different auxiliary spaces. 

Moreover, jz/T*) is defined as |z/z) in Eq. (I14L but with p* 
instead of p. Similarly, rt^{p*) is defined in a way similar 
to Eq. ([7D; 

n'^ip*) := + T+a- + T_a+ 


where the operators Ta are defined with p* instead of p. 
Moreover, they commute with the Sa since they act on 


different auxiliary spaces. 

Next define 

Ha = tr[cr“H(p)H'^(p*)], a€{0,x,y,z}. (15) 

Explicitly we have 

Bo = 2S^T, + S+T+ + S-T_ (16) 

H, = (H_ - S+)T, + S,{T_ - T+) (17) 

By = i[S,{T_ + T+) - {S. + S+)T,] (18) 

B, = S+T+ - S-T_ . (19) 

The spin current may then be written as 




and that — Sz commutes with Bq. This reflects the 
translational symmetry of the Ji in the steady-state. 
Hence, making use of Eq. dia, we arrive at 


27 Z{N - I) 
72 -p h2 Z{N) 


( 21 ) 


This is the required formula for the steady-state magne¬ 
tization flux. In the present model J is a function of TV, 
h, 7 and 9 only. Eq. (EID must be computed for each TV. 
Even though this may be done exactly, the formulas be¬ 
come extremely cumbersome for large sizes. On the other 
hand, computing J numerically is now a trivial task. 

Also of interest is the magnon density irii) = (1-1- 
(crf))/2. A calculation similar to the above leads to 


K) 


Z{N) 


( 22 ) 


IV. RESULTS 

We now discuss the behavior of J as a function of TV, 
h, 7 and 6. The focus will be on the case / = 1, for which 
the MPS solution is valid. Notwithstanding, a few words 
will be given about the case / < 1. 

We begin with 9 = 0 and h = 0. The spin current as 
a function of TV and 7 is presented in Fig. [ 2 ] In order 
to interpret these results, recall that magnons are con¬ 
stantly being pumped at the left source, which then prop¬ 
agate through the lattice and are eventually collected in 
the right drain. The spin current is then simply pro¬ 
portional to the number of magnons being collected at 
the right drain. This number depends on two things: (i) 
the number of magnons being injected per unit time in 
the left source, which is proportional to 7 and (ii) the 
magnon scattering events during the trip to the right 
drain. In standard electrical conduction (e.g. in Drude’s 
model), the electrons scatter with lattice imperfections 
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FIG. 2. (Color Online) Spin current J for f = I, 9 = 0 and 
h = 0. (a) J vs. 7 for different sizes N. (b) J vs. N for 
different values of 7 . The dotted black line has slope -2. 


or phonons. Since the number of scattering agents scales 
proportionally to N we then have a diffusive current 
J ^ 1/N. In our case the magnons do not scatter with 
lattice imperfections. They either travel through unim¬ 
peded or they participate in 4-magnon scattering events 
(where 2 magnons scatter producing two new magnons 
in the process^^). When 7 is sufficiently small the den¬ 
sity of magnons in the chain is very small, thus making 
these events very rare. In this case J will increase with 7 
and will also be independent of N; i.e., ballistic. This is 
clearly observed in Fig.j^Ka), where we see that the curves 
for different N overlap when 7 is small. Conversely, in 
the high 7 limit the number of magnons, and hence the 
number of scattering events, will be significant. In this 
regime it is found^^ that J is sub-diffusive, behaving as 
J ^ l/N'^. The reason for this is that by doubling the size 
of the chain, we quadruple the number of four-magnon 
scattering events. As shown in Ref. the transition 
between the ballistic and sub-diffusive regimes occurs at 


7 


1 

N 


(23) 


A clear example of this transition is seen in the curve for 
7 = 10“^ in Fig.[2](b), where the regime changes abruptly 
from J ^ I to J N~^ exactly at N = 100. 

Next we discuss the behavior for non-vanishing bound¬ 
ary fields, h ^ 0, still keeping 9 = 0. In Fig. [SKa) we 
present N'^J vs. h for 7 = 1 (sub-diffusive; high magnon 
density). As can be seen, even for moderately small sizes, 
the curves start to scale very well according to J ^ 1 /A^^. 
In this scaling region we have found that the current is 
very well described by 


which is illustrated by the dashed line in Fig. |3l(a). Note 
also that J is asymmetric with respect to h' i.e., the spin 
current is rectifiedJ^ 

The changes which occur as we reduce 7 below 7 * 
are illustrated in Fig. EKb), where we plot J /7 vs. h 
for N = 100 and different values of 7 . As can be seen, 
there is a drastic behavioral transition from a bell-shaped 
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FIG. 3. (Golor Online) Spin current J as a function of the 
boundary fields h with / = 1 and 9 = 0. (a) N^J vs. h for 
7 = 1 and different values of N. (b) J /7 vs. h for N = 100 
and different values of 7 around 7 * = 1/A = 0.01. (c) and 
(d) J /7 vs. h for 7 = 10~® and different values of A. The 
dashed lines in (d) correspond to Eq. (I24II . 


structure at 7 > 7 * to a plateau at 7 < 7 *. This plateau 
is illustrated in more detail in Figs. ISKc) and (d) for 
7 = 10“^ and different sizes. As can be seen, the plateau 
region is asymmetric with respect to h and independent 
of size. It corresponds to the ballistic behavior of the spin 
current. As the field is increased, however, one eventu¬ 
ally observes an abrupt transition to a much lower spin 
current. For positive fields the transition is continuous 
whereas for negative fields it is discontinuous (strictly 
speaking, it is only discontinuous in the thermodynamic 
limit). The critical field where the plateau transition oc¬ 
curs is found from the simulations to be h* ~ —5/A. We 
also call attention to the fact that outside the plateau re¬ 
gion, J is again well described by Eq. (HI, as illustrated 
by the dashed lines in Fig. (H^d). This indicates that for 
large fields the behavior is again sub-diffusive. 

The results presented so far were obtained from the 
exact MPA steady state which is valid only at / = 1 (zero 
temperature). However, the rich behavior of the current 
observed for / = 1 also survives at finite temperatures; 
i.e., for / < 1. This can be seen in Fig.|4]where we report 
the current J vs. / as obtained from the exact numerical 
diagonalizationii of Eq. m for A = 6 . The current as 
seen from the numerics shows basically the same features 
as in the MPA case: a bell-shaped behavior at high 7 and 
a sharp plateau at low 7 (for this small size the plateau 
is not yet completely formed). In Eig. 0] we also plot 
the MPS solution when / = I to illustrate the perfect 
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FIG. 4. (Color Online) Spin current J vs. h when / < 1, 
computed using the exact diagonalization of Eq. ([2| for N = 
6 . (a) 7 = 10~®. (b) 7 = 1. The dotted black lines correspond 
to the MPS solution when / = 1. 
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FIG. 6 . (Color Online) J vs. h for N = 500 and different 
values ol 9 (as defined in Fig. 0. (a) 7 = 10 ^ and (b) 7 = 1 . 


FIG. 5. (Color Online) Small size effects in the spin current, 
(a) J /7 vs. Nh for 7 = 10~® and different values of N. (b) 
J /7 vs. Nh for = 15 and different values of 7 . 


agreement between both methods. 

The gradual formation of the plateau as the size of the 
system increases in illustrated in Fig. (SKa). In Fig. [S{b) 
we show the changes which occur as one changes 7 when 
N = 15. As can be seen in both images and in Fig. He), 
when N is small the current presents a series of irregular 
and sharp resonances when /i < 0 , at positions which 
vary with N (such peaks have been observed recently in 
Ref. [33) . It is important to note, however, that these 
peaks only appear for 7 < 1 /iV^ and therefore become 
vanishingly small for any moderately large size. This 
can be seen, for instance, by comparing the curves with 
7 = 10“^ and 7 = 10“^ in Fig.jSKb). Both are practically 
identical, except for the peaks, which are only present 
when 7 = 10“^. Note also that it follows from Eq. m 
that J is bounded so that these cannot be delta peaks. 

We consider now the case with a general twisting angle 
0 € [0, tt]. Fig. [HKa) shows J vs. h for different values of 
9 with fixed size N = 500 and 7 = 10“^. As expected, 
J — ^ 0 as 0 — > TT. However, and remarkably, even for 
values of 9 close to the undriven situation 0 = tt, one 
still observes high values of J for negative values of ft., 
in a plateau region that shrinks as ft ^ tt. Thus, by 
monitoring the twisting angle ft, one can fine-tune the 
high current plateau width. For completeness, we also 
show the behavior for large 7 in Fig. Hb). 



Lattice site i 



Lattice site i 


FIG. 7. (Golor Online) Magnon density profile (rii) = (1 -|- 
(crf))/2 for different values of ft, with N = 500, 7 = 10~® and 
ft = 0. (a) ft < 0 near the plateau transition [cf. Fig. [3jd)]. 
(b) ft > 0 . 


V. DISCUSSION AND CONCLUSIONS 


The remarkable and sharp transitions observed in the 
spin current, from ballistic (inside) to sub-diffusive (out¬ 
side the plateau), as the magnitude |ft| of the boundary 
fields is increased, suggest that sufficiently high fields act 
as scattering barriers, impeding magnons to flow through 
the system, from source to drain. This can also be seen by 
looking at the magnon density profile (n^) = (1 -I- ((Tf))/2 
plotted in Fig. [7] for N = 500, 7 = 10“^ and ft = 0. The 
red (solid) curve in Fig. Ha) corresponds to the profile 
in the plateau (ballistic) region of Fig.^d). In this case, 
the distribution is flat with (n^) ~ 1 / 2 , characteristic of 
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a maximal current state. On the other hand, outside the 
plateau the profile is sine-shaped, characteristic of the 
sub-diffusive regime.— The transition between the two 
profiles is discontinuous for h < 0 [Fig. |71[a)] and contin¬ 
uous for h > 0 [Fig. [2l[b)]. Hence, we conclude that the 
density of magnons inside the chain may also be adjusted 
by changing the boundary field h. Chumak et. al^ used 
a similar idea to construct their magnonic logic gate. But 
in their case an additional source of magnons was respon¬ 
sible for changing the magnon current and the magnon 
density. Consequentially, the transition between the on 
and off states was in their case quite smooth. Here we 
see an extremely abrupt transition, thus being potentially 
more suited for a logic gate. 

In what concerns an experimental realization of the 
present idea, it is important to note that even though we 
studied a very specific situation, the underlying physi¬ 
cal principles of our results are very general, being based 
only on the entrapment of magnons by magnetic fields. 
Hence, similar results should be obtained in different 
field configurations which maintain the same principles. 
Most magnonic circuits are constructed using Yttrium 
iron garnet (YIG),— which is well described by the 
Heisenberg model, albeit with a different spin value. The 
Lindblad generators then represent microstrip antennas 
which are used to generate and collect magnons.— 
Even though the Lindblad dissipators have been exten¬ 
sively used in the past to study open quantum systems, 
we are unaware of any papers mentioning this specific 


application of them as describing the injection and col¬ 
lection of magnons. 

The energy and time units of the problem are set by 
the constant which should appear in the first term of 
Eq. ([T]), but which we have throughout set as unity. Ac¬ 
cording to Ref. [s^ J ~ 10“^^ J. The pumping rate 7 
(measured in magnons per second) should operate below 
the critical value 7 * which, in the correct units, reads 
7 * = J/hN ^ 10^^/A^ Hz. This gives the optimal value 
of 7 below which the flux should be ballistic. Letting h = 
/isH, where is the Bohr magnetron, we find that the 
critical magnetic field B* where the plateau transition 
occurs is, in correct units, |i?*|(T) ~ J/fisN ~ lO/N. 
Hence, for any reasonable values of very small mag¬ 
netic fields may suffice to induce the plateau transition. 

In summary we have studied the quantum Heisenberg 
chain driven by two Lindblad baths and subject to two 
magnetic fields acting on each boundary. An exact solu¬ 
tion was given in terms of matrix product states which 
enables one to calculate local observables for any chain 
size. The system is seen to undergo a discontinuous tran¬ 
sition from ballistic to sub-diffusive spin current as a 
function of the field intensity. Thus, the system may 
function as an extremely sensitive magnonic logic gate 
using the boundary helds as the base. 
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The Heisenberg Hamiltonian prohibits 3 magnon events 
since it conserves the number of excitations.— Events in¬ 
volving more than 4 magnons are in principle possible, but 
much less likely. In real systems 3-magnon events do exist 
M)'a305n^I^lfi:):ice of more complex interactions (e.g. dipo¬ 
lar). However, as shown in Ref. 4-magnon events are 

the most relevant type, at least for YIG. 
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